clear;
G=[1,0,0;0,3,0;0,0,3];
x0=[7;8;9];
f0=1/2*x0'*G*x0;
g0=G*x0;
d0=-g0;
alpha0=-g0'*d0/(d0'*G*d0);
x1=x0+alpha0*d0;
f1=1/2*x1'*G*x1;
plot3(x0(1),x0(2),x0(3),'ro');
hold on;
[X Y Z]=ellipsoid(0,0,0,sqrt(2*f1/G(1,1)),sqrt(2*f1/G(2,2)),sqrt(2*f1/G(3,3)),30);
surf(X, Y, Z,'FaceAlpha',0.1)
plot3(x1(1),x1(2),x1(3),'ro');
plot3([x0(1), x1(1)],[x0(2), x1(2)],[x0(3), x1(3)],'r','LineWidth',2);
g1=G*x1;
beta1=g1'*G*d0/(d0'*G*d0);
d1=-g1+beta1*d0;
alpha1=-g1'*d1/(d1'*G*d1);
x2=x1+alpha1*d1;
plot3(x2(1),x2(2),x2(3),'r*');
f2=1/2*x2'*G*x2;
[X Y Z]=ellipsoid(0,0,0,sqrt(2*f2/G(1,1)),sqrt(2*f2/G(2,2)),sqrt(2*f2/G(3,3)),15);
surf(X, Y, Z,'FaceAlpha',0.3)
plot3([x1(1), x2(1)],[x1(2), x2(2)],[x1(3), x2(3)],'r','LineWidth',2);
g2=G*x2;
beta2=g2'*G*d1/(d1'*G*d1);
d2=-g2+beta2*d1;
alpha2=-g2'*d2/(d2'*G*d2);
x3=x2+alpha2*d2;
plot3(x3(1),x3(2),x3(3),'r*');
plot3([x2(1), x3(1)],[x2(2), x3(2)],[x2(3), x3(3)],'r','LineWidth',2);
axis equal;
